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ABSTRACT 

The process of gravitational scattering of planetesimals by a massive protoplane- 
tary embryo is explored theoretically. We propose a method to describe the evolu- 
tion of the disk surface density, eccentricity and inclination caused by the embryo- 
planetesimal interaction. It relies on the analytical treatment of the scattering in two 
extreme regimes of the planetesimal epicyclic velocities: shear-dominated (dynamically 
"cold") and dispersion-dominated (dynamically "hot"). In the former, planetesimal 
scattering can be treated as a deterministic process. In the latter, scattering is mostly 
weak because of the large relative velocities of interacting bodies. This allows one to 
use the Fokker-Planck approximation and the two-body approximation to explore the 
disk evolution. We compare the results obtained by this method with the outcomes 
of the direct numerical integrations of planetesimal orbits and they agree quite well. 
In the intermediate velocity regime an approximate treatment of the disk evolution is 
proposed based on interpolation between the two extreme regimes. We also calculate 
the rate of embryo's mass growth in an inhomogeneous planetesimal disk and demon- 
strate that it is in agreement with both the simulations and earlier calculations. Finally 
we discuss the question of the direction of the embryo-planetesimal interaction in the 
dispersion-dominated regime and demonstrate that it is repulsive. This means that the 
embryo always forms a gap in the disk around it, which is in contrast with the results of 
other authors. The machinery developed here will be applied to realistic protoplanetary 
systems in future papers. 

Subject headings: planets and satellites: general — solar system: formation — (stars:) 
planetary systems 



1. Introduction. 

This paper continues the line of investigation started in our previous work (Rafikov 2001, 
2002a; hereafter Papers I and II) which was devoted to the treatment of planetesimal-planetesimal 
gravitational interactions. Here we consider the interaction between the growing protoplanetary 
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embryo and the planetesimal disk. By embryo we imply in the present context a single body 
with mass M e much larger than the masses of individual planetesimals m. There are several 
reasons for studying this important problem separately from the mutual gravitational scattering of 
planetesimals. 

First, planetesimal-planetesimal encounters in realistic protoplanetary disks usually occur in 
the dispersion-dominated regime, which applies when the relative approach velocity of two particles 
is bigger than the differential shear in the disk across the Hill (or tidal) radius. The Hill radius is 
defined as 



where a is a value of the semimajor axis at which the interaction takes place, and mi, rri2, and M c 
are the masses of interacting planetesimals and of the central star. 

However, in the same protoplanetary disk gravitational interaction between the embryo and 
planetesimals could be in the opposite velocity regime — shear-dominated — when the planetesimal 
random motion is small compared to the shear across the Hill radius, simply because the embryo 
mass and therefore the Hill radius is much larger. Indeed, in the case of embryo-planetesimal 
interactions Hill radius Rh = a e (M e /M c ) 1 / 3 3> rn (here o e is a semimajor axis of the embryo) 
as a consequence of M e » mj. Thus, reduced (normalized in Hill coordinates) values of random 
velocities in the embryo-planetesimal case are smaller by a factor [(mi + m2)/M e ] 1 / 3 <C 1 than 
those corresponding to the planetesimal-planetesimal interactions. 

Of course, it might be that the planetesimal disk has already been so excited dynamically that 
even embryo-planetesimal encounters are in the dispersion-dominated regime. Thus, we consider 
here both regimes of the embryo-planetesimal interaction — shear- and dispersion-dominated. 

Second, embryo-planetesimal interactions are complicated by the presence of a special type of 
orbits in a 3-body problem — the so-called horseshoe (or librating) orbits (Henon & Petit 1986; 
Murray & Dermott 1999). Planetesimals on these orbits do not perform the usual circulating mo- 
tion which is characteristic of passing orbits (the most important case for planetesimal-planetesimal 
scattering) but a librating one. This horseshoe motion can only occur when the difference in semi- 
major axes of interacting bodies is smaller than their Hill radius. For planetesimal-planetesimal 
interactions th is negligible compared to the scale of surface density variations or the radial epicyclic 
excursion of an individual planetesimal. Thus horseshoe orbits are unimportant in this case. How- 
ever, the Hill radius of the embryo-planetesimal interaction Rh can be comparable to the length 
scale of the disk inhomogeneities caused by the embryo. Thus the phenomenon of horseshoe motion 
can be crucial for the planetesimal dynamics near the embryo (see below §3.1). 

Third, as we have already mentioned in Paper II, planetesimal-planetesimal scattering is de- 
scribed in terms of the disk properties averaged over some region of the disk, which diminishes the 
importance of the details of spatial distributions of disk properties. On the contrary, in the case of 




(1) 



-3- 



the embryo-planetesimal interaction we are interested in details of the spatial behavior of various 
quantities characterizing the state of the disk and they are the primary goal of our present study. 

All these complications preclude the direct application of the results obtained in Paper II to the 
present consideration. However the general analytical approach to the treatment of planetesimal 
disk evolution developed there remains valid and we will employ it in this paper. 

Numerical orbit integrations (Tanaka & Ida 1996, 1997) and N-body simulations (Ida & Makino 
1993) provide an alternative and important route of studying embryo-planetesimal interactions. 
Their drawback is their intrinsically low speed and inability to treat large number of planetesimals. 
However, since the physics is incorporated in them on a very basic level with the minimum of 
additional assumptions they can provide us with robust predictions. To use this advantage of 
numerical methods and to avoid their handicaps we employ the following methodology: we provide 
a self-consistent analytical description of the embryo-planetesimal interaction in different velocity 
regimes. To check this description and to verify the validity of the simplifying assumptions utilized 
in its development we have used numerical orbit integrations performed for several sets of typical 
planetesimal disk parameters. After we make sure that our theory works well for these sets of 
parameters we can use it for others as well and be confident of its reliability. What we gain by this 
approach is the speed of computation and ability to explore the whole space of important physical 
variables. 

The condition on the embryo's mass, M e S> m, has important dynamical implications. In many 
applications it justifies the neglect of the embryo's recoil resulting from planetesimal scattering. 
Also, dynamical friction between the embryo and planetesimals will tend to produce random energy 
equipartition (Stewart & Wetherill 1988; Wetherill & Stewart 1989) which means that embryo's 
eccentricity and inclination are most likely to be negligibly small. Thus, we will assume in this paper 
that the embryo moves on a fixed circular orbit and its eccentricity and inclination are zero. We will 
also consider the embryo to be isolated from the gravitational effects of other massive bodies which 
may be growing nearby, an assumption which can easily be abandoned in future work. Throughout 
the paper we neglect the presence of any resonant effects. This is justifiable if frequent planetesimal- 
planetesimal encounters can destroy any commensurabilities with the embryo's rotation period. The 
embryo's recoil in the course of planetesimal scattering and distant embryo-embryo interactions 
would also help to do that. We leave for the future the clarification of conditions necessary for 
employing this simplification. 

It will be convenient to use the embryo's Hill radius Rh as a unit of length in our study. We 
introduce the Hill orbital elements of a planetesimal evaluated at large azimuthal distance from 
the embryo — the difference in semimajor axes H, eccentricity e, and inclination i relative to the 
embryo's orbit at semimajor axis a e : 



H = 



a — a e 




(2) 



Rh 



and use them as planetesimal coordinates. The distribution function of the orbital elements e and 
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i will be assumed to have a Rayleigh form with dispersions a e and af. 

(3) 



ip(e, ijaeai = ~ 2 ~2 ex P 



Tal ~ 2a? 



We will also be using the dimensionless surface number density of guiding centers N(H) to char- 
acterize the planetesimal spatial distribution. 

Although we focus on a single embryo, for many applications we can treat the embryo using a 
continuous form of the evolution equations. Fro example, we may assume that the discrete surface 
number density of the embryo is given by 

N em (H) = -U(m - M e )5(H). (4) 

Since planetesimal-planetesimal interactions are not important here it will be enough to consider a 
single-mass planetesimal population. 

The three-body interaction in the Hill approximation preserves a certain combination of relative 
orbital elements of interacting bodies called the Jacobi constant (Goldreich & Tremaine 1980; Henon 
& Petit 1986): 

J = e 2 + i 2 -^F 2 + 20 e , (5) 

where (fr e is the gravitational potential of the embryo, which can be neglected far from the embryo. 
For embryo-planetesimal scattering one can introduce the concept of integrated Jacobi constant of 
the whole planetesimal population: 



J 



tot 



oo 

/ 



2N{H)a 2 e + 2N{H)af - jN(H)H 2 



dH. (6) 



This quantity should be conserved because (1) each individual planetesimal scattering off the em- 
bryo conserves the Jacobi constant of the relative motion, and (2) embryo's random motion is 
negligible, which means that relative eccentricity, inclination, and difference in semimajor axes are 
determined by planetesimal orbital parameters only. We will use the conservation of this quantity 
and of the total number of planetesimals (we neglect their coagulation at this point) 

oo 

N tot = f N ( H ^ dH ( 7 ) 



as checks of our evolution equations. 

The orbit integrations that we use are performed by solving Hill equations numerically. We 
have integrated the evolution of the system [equations of osculating orbital elements evolution 
(11) of Paper II] using fourth order Runge-Kutta integrator (Press et al. 1988). Unlike similar 
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calculations of Tanaka & Ida (1996, 1997) our orbit integrations do not employ additional analytical 
simplifications to avoid possible biases. In a typical integration the Jacobi constant is conserved 
with fractional accuracy 1CT 8 — 1CT 12 . The results of these orbit integrations and their comparisons 
with theoretical predictions will be presented in the following sections. 

We devote §2 to studying the shear-dominated case and §3 to exploring the dispersion-dominated 
case. The velocity regime intermediate between them is addressed in Appendix B. We discuss some 
general features of the embryo-planetesimal interaction in §5. Some auxiliary results are presented 
in appendices: Appendix A contains the derivation of the probability distribution of scattered semi- 
major axes in the dispersion-dominated regime, while in Appendix C we calculate the embryo's 
accretion rate in different velocity regimes. 



2. Scattering by the embryo in the shear-dominated regime. 

We consider first the embryo-planetesimal interaction in the shear-dominated regime. In Pa- 
per II we have derived a master equation (30) for the evolution of the planetesimal distribution 
function. We will use this general equation and equation (33) of Paper II as a basis for our further 
developments. The conditions for the shear-dominated regime to be realized are 

d 2 e < 1, af < 1. (8) 

Some important simplifications can then be made. 

First, scattering in this regime is deterministic in the sense that the outcome of an interaction 
between two particles depends only on their difference in semimajor axes before the collision Hq 
and not on their relative random motion (which is naturally absent in this case). This means that 
the change of the reduced semimajor axis difference Ahgc 1 is a single- valued function of only Hq in 
this regime. 

Second, inclination is hardly excited at all in the course of an encounter. A considerable change 
of inclination requires a substantial force acting in the vertical direction. However, the dynamically 
cold disk is very thin and it is easy to see that the gravitational force between the interacting 
particles is directed almost horizontally. Thus, noticeable inclination growth can occur only for 
particles experiencing very close encounters. From qualitative considerations one would expect 
that change of the inclination vector i (see Petit & Henon 1986; Ida 1990; Paper II) is given by 

A(i sc ) 2 = i2 5l (F )«l, (9) 

where io is the initial value of the reduced vector inclination and g\{Ho) ~ 1 is some function which 
can be easily computed numerically. For our purposes we neglect the inclination growth due to the 
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gravitational stirring in the shear-dominated regime completely and set Ai sc = 0. We will however 
keep the terms describing the transport of vertical energy in the disk (see below). 

The absence of heating in the vertical direction naturally leads to another simplification. Since 
we can neglect the change of inclination in the encounter the change of eccentricity becomes directly 
related to the change of semimajor axis difference due to the conservation of Jacobi constant (5). 
Then we can write that the change of the vector eccentricity e is 

A(e sc ) 2 = |a(/&) = \ [{Ah sc f + 2H Ah sc ] . (10) 

Note that A(e sc ) 2 ~ 1 if Hq ~ 1. It can also be easily shown that 

eAe sc ~e 2 «l, i • Ai sc ~ i 2 < 1 (11) 

in the shear-dominated regime; the corresponding scattering coefficients in evolution equations will 
therefore be neglected in this paper. 

Thus, the probability P r of scattering from Ho,eQ,iQ to H, e, i can be written in the shear- 
dominated regime as 

P r = 5[AH - Ah sc (H )} S[Ae - Ae sc (H )] S(Ai), (12) 

where S denotes the Dirac delta function. 

The computational challenge of the shear-dominated regime is that strong scattering is possible 
for Hq ~ 1, i.e. AH in the course of an encounter can be quite substantial. Thus we cannot use the 
Fokker-Planck formalism for the shear-dominated scattering and one has to deal with the different 
moments of the master evolution equation (30) of Paper II in their general form. However, this 
does not pose an insurmountable problem because the probability distribution function of the 
shear-dominated scattering is a single- variable function only. An analytical fit to this function was 
calculated by Petit & Henon (1986, 1987b) using results of numerical orbit integrations. This fit 
automatically takes horseshoe motion into account so that we do not have to worry about the 
complications associated with this type of orbit in the shear-dominated regime: one can see from 
the expression for Ah sc (Ho) (Petit & Henon 1987b) that Ah sc (Ho) — ► —2Hq in the shear-dominated 
regime as Hq — ► 0. 

Using the deterministic form of P r [substituting expression (12) into equation (30) of Paper II] 
and the embryo's surface density in the form (4) we can take different moments of e 2 and i 2 . As a 
result we find the following set of equations describing the evolution of the planetesimal population 
due to the scattering by the embryo in the shear-dominated regime: 



dN(H) \A\fi 



1/3 
e 



dt TT 



| [2N(ff)o*{H)] = - 



N(H)\H\- J P(H ,H)N(H )\Ho\dH 

— oo 

2N{H)d 2 e {H)\H\ 



(13) 
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j P(H ,H)N(H ) (2a 2 e (H ) + A(e sc ) 2 (H )) \H \dH 



(14) 



d [2N(H)d?{H)]=-\ Alf ''"' 



dt 



TT 



2N(H)a 2 (H)\H\ 



P(H ,H)2N(H )a 2 (H )\H \dH 



(15) 



with n e = (Me/Mc) 1 / 3 . Here A is the Oort's constant characterizing the differential rotation of 
the disk (Binney & Tremaine 1987), A = -(3/4)0 in a Keplerian disk (Q, = ^/GM c ja\ is the disk 
rotation frequency) , Hq is an integration variable having the meaning of the initial difference of the 
semimajor axes of interacting bodies, A(e sc ) 2 is defined in equation (10), and 



P(H , H) = 5 H-Hq- Ah sc (H ) 



(16) 



Terms in the r.h.s. of equation (15) are ~ a, <C 1, i.e. are of the same order as A(i sc ) 2 which 
we have agreed to neglect. This inconsistency stems from the fact that we want our evolution 
equations to preserve the integrated Jacobi constant (6). For this reason we keep the terms in the 
r.h.s. of (15) (which are essentially the transport terms) but neglect terms like A(i sc ) 2 if A(e sc ) 2 
is assumed to be given by (10). Then dJ tot /dt = which can be verified using equations (10) and 
(13)-(15). Also, one can easily check that equation (13) conserves the total number of planetesimals 
in the disk. 

One technical issue merits mentioning at this point. In the shear-dominated regime the semi- 
major axis difference after the encounter H is a single- valued function of Hq. However, the inverse 
function Hq(H) is multivalued (Petit & Henon 1987b). For this reason integrals of some quantity 
F(H ) [e.g. F = N(H )\H \] over P(H ,H)dH in (13)-(15) result in 

r ~ i -i 

dAh 

P(H ,H)F(H )dH = YF(H ok ) ^ 



/ 



k 



1 + 



dH 



where i^ofc is the k-th root of the equation 2 

H = H + Ah sc (H ). 



(17) 



(18) 



The system (13)-(15) forms a closed set of equations needed to describe the disk evolution 
caused by embryo-planetesimal scattering in the shear-dominated regime. Equation (13) has already 
been derived by a different method in Paper I, and now we have extended that analysis by taking 
random velocity evolution of the disk into account. Following Paper I we will be using AJi sc (Hq) 
in the analytical form suggested by Petit & Henon (1987b; see also Appendix B of Paper I). 



2 This form arises because of the integration rule of the Dirac 5-function: f p(x)8[q(x)] = ^2p{xj)/q'(xj), where 

-oo i 

Xj is the j-th root of equation q(x) = 0. 
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3. Scattering by the embryo in the dispersion-dominated regime 

When the embryo-planetesimal scattering is in the dispersion-dominated regime the same 
simplifications as in the case of planetesimal-planetesimal scattering can be made: scattering is weak 
and this warrants the use of both the two-body approximation and the Fokker-Planck expansion. In 
this regime we can use many of the results obtained in Paper II. There are however some additional 
complications; one of them is related to the aforementioned existence of horseshoe orbits. To discuss 
this issue we need to know which conditions should be fulfilled for horseshoe motion to take place 
in the dispersion-dominated velocity regime. 

3.1. Horseshoe orbits. 

Horseshoe motion arises when the separation of the semimajor axes of interacting bodies is 
small. In this case their relative motion due to the shear in the disk is very slow; for this reason 
even the weak gravitational force between the bodies at large azimuthal distance can lead to a con- 
siderable change of their angular momenta over a long time interval; the relative motion of particles 
reverses, they turn around and move on almost closed trajectories until the next conjunction (see 
Figure la). In the Solar System the best example of such motion is given by Saturnian satellites 
Janus & Epimetheus (Murray & Dermott 1999). 

In contrast, planetesimals on pasing orbits pass each other after the interaction, usually without 
reversing their motion (see Figure lb). In the shear-dominated regime these orbits typically have 
initial separations of semimajor axes larger than the Hill radius. However in dispersion-dominated 
encounters the question of the spatial separation of horseshoe and passing orbits becomes more 
subtle. 

Namouni (1999) has considered dispersion-dominated scattering in the planar case (i = 0) 
and suggested that for a given relative eccentricity e> 1 the boundary of the horseshoe region is 
located at H = hh s ~ e -1 / 2 [see also Ida & Nakazawa (1989)]. This means that the more "energetic" 
the planetesimals are, the smaller their semimajor axis separation has to be for horseshoe motion 
to appear. This tendency is illustrated in Figure 1: orbits with the same value of H = 0.7 
exhibit different behavior depending on their random velocities — in the shear-dominated regime 
(e = i = 0) horseshoe motion takes place while in the dispersion-dominated (e = i = 5) the orbit 
is passing; however if one decreases H this orbit also becomes horseshoe (which is shown by a 
trajectory with H = 0.3, e = i = 5 in Figure la). It is easy to understand the reason for this 
dependence on the planetesimal random velocity: the higher the eccentricity the more time the 
approaching planetesimal spends far from the scatterer because it moves on a very large epicycle. 
As a result the mutual force is weak, the angular momentum exchange is small, and the planetesimal 
passes the scatterer instead of turning around. In other words, the higher e the less noticeable the 
presence of the scatterer is for the incoming planetesimal. 
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Fig. 1. — Representative trajectories of bodies on horseshoe (a) and passing (6) orbits for different 
values of initial orbital elements. The scattering body is at the origin and the coordinates are in 
units of the Hill radius. 

In the nonplanar case, on the basis of similar arguments, we suggest that the horseshoe bound- 
ary condition in the high-velocity regime should be replaced with 

/ k \ 1/4 

A "»=(iiT?) ' (19) 

To check this prediction and fix the proportionality constant k, we have performed a set of orbit 
integrations in both the shear- and dispersion-dominated regimes. We have separated the outcomes 
of encounters into horseshoe orbits (when H was changing sign as a result of an encounter) and 
passing ones (when the sign of H remained unchanged). The results are presented in Figure 2. One 
can see a rather clean separation between the horseshoe orbits (red dots) and passing ones (blue 
dots). There are some red dots which appear in the region mostly occupied by the passing orbits, 
but they originate from large-angle scattering during close encounters and are not librating orbits 
as normal horseshoes are. 



0.01 



0.1 

H 



1 



Fig. 2. — Horseshoe and passing orbits as a function of the orbital parameters H and e 2 + i 2 . Red 
ones correspond to horseshoe orbits, blue ones to passing orbits. The solid line is the analytical 
expression (20) for the boundary. The dashed line shows the power law asymptote of this relation 
in the strongly dispersion-dominated regime. 

The separation condition (19) fits the boundary between the two types of orbits very well 
when e 2 + i 2 S> 1 if we take k — 8 (dashed curve in Figure 2). However, as we move away from 
the strongly dispersion-dominated regime some deviations from (19) appear. This is not surprising 
because small e and i correspond to the shear-dominated regime, and then the horseshoe region 
has a well-defined boundary at hhs ~ 1- We have found that the shape of the horseshoe-passing 
boundary can be fit rather accurately by the following condition 



which is shown by the solid curve on Figure 2 3 . 



3 Our value of b would predict hhs in the shear-dominated regime slightly different from that suggested by Petit & 
Henon (1987b): ~ 1.4 instead of ~ 1.2. However, such a small difference is unlikely to be important for our purposes. 
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1.4 



(20) 
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Using distribution function of eccentricities and inclinations (3) we can calculate the fraction 
of passing orbits p paS s- 



Ppass (-^0 




+ 



exp 



2uf 



— exp 



H < b, 
H>b. 



(21) 



One can see that p pa ss{H) changes rapidly from to 1 between the regions of horseshoe and 
passing orbits (this is caused by the exponential dependence of p pass on Rh s and the strong power 
law dependence of Rh s on H). This allows us to neglect the transition region in the dispersion- 
dominated regime and assume these different types of orbits to occupy different regions of the 
physical space. Thus we take passing orbits to be restricted to the region ^> hhs where hh s is 
defined such that 

Ppassihs) = 1/2, (22) 

and horseshoe orbits to satisfy \H\ < h^g. Using these conditions we can now consider passing 
orbits separately from the horseshoe ones. 



3.2. Scattering on horseshoe orbits. 

The relative motion of interacting planetesimals in the horseshoe regime can be split into the 
slow shear motion of the guiding centers and the fast epicyclic motion. When such a separation 
is legitimate some quantities associated with the fast motion called adiabatic invariants should be 
conserved (Landau & Lifshitz 1989). It was shown by Henon & Petit (1986) and Hasegawa h 
Nakazawa (1990) that the absolute value of the relative semimajor axis difference, relative eccen- 
tricity and relative inclination are all separately conserved quantities in the course of a horseshoe 
encounter, and that they are the adiabatic invariants of this type of motion. Thus the quantities 
//o,eo,io before the encounter are related to their values after, H, e, i, as 

H = -H , e = e , i = i (23) 

It follows from this conservation of adiabatic invariants that the effect of the embryo-planetesimal 
encounters in the horseshoe region is to exchange planetesimals at symmetric orbits (H and — H). 
Interacting bodies just librate between successive close approaches when they reverse the direction 
of their motion. 

Because the embryo is much more massive than the planetesimals and its e and i are zero, 
the relative motion in the embryo-planetesimal system is equivalent to the motion of planetesimals. 
The scattering probability for the horseshoe motion i\ s looks like i\ s = S(AH + 2Ho)5(Ae)5(Ai) . 
Using this result we can write down the evolution equation of some quantity F by analogy with 
equations (13)-(15) in a very simple form: 

^J*f m[Fi - H) - Fm . (24) 
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In our case F can be the planetesimal surface density N, horizontal epicyclic energy Na 2 or vertical 
epicyclic energy No 2 



r 2 e 



3.3. Scattering on passing orbits. 

To describe the embryo-planetesimal interaction on passing orbits in the dispersion-dominated 
regime we will use the results for the scattering in this regime derived in Paper II, namely equations 
(49)-(51) and (52)-(55). We substitute embryo's surface number density in the form (4) for the 
surface density of approaching bodies N2, assume that m/M e — > and that a e 2 = <5"i2 = 0. As a 
result we obtain the following system 



~(|HK^) + ^(|H|<(M)V) 



dN _ \A\t/ e /3 

Ot 7T 

I [2N(H)a 2 e (H)] = 

\H\(A(e 2 ))N - A (|F|((e 2 + 2e • Ae)A~h)N) + \ A_ (\H\(e 2 (A~hf)N 



(25) 



(26) 



and an equation for the inclination evolution analogous to (26). These formulae are only valid 
in the region of space restricted by the condition (22). Analytical expressions for the scattering 
coefficients (Ah), ((Ah) 2 , (A(e 2 )), ((e 2 + 2e • AS) Ah), and (e 2 (Ah) 2 ) can be found in Paper II. The 
proper boundary conditions for this system will be derived in §3.3.1. 

At this point we should address a subtle issue which was not important for the planetesimal- 
planetesimal scattering but becomes nontrivial for the embryo-planetesimal interaction. It is related 
to the fact that planetesimals at different initial separations are driven past the embryo by the 
differential shear at different rates. Those which initially had \H\ 3> 1 quickly approach the 
embryo, experience scattering, and quickly depart. On the approach and departure stages their 
orbital parameters do not have time to change except when they are close to the embryo. Thus, 
we can assume that values of H,e,i far from the embryo are the same as at the point of closest 
approach if initially \H\ 3> 1. 

However, planetesimals which start their motion not far from the horseshoe orbits boundary 
are not moving very fast (because this boundary is located at hhs < Rh)- As a result, for these 
planetesimals the exchange of the angular momentum before and after the close encounter can be 
important, just as in the case of horseshoe orbits, where such an exchange essentially prevents close 
encounters from happening. Consequently, the asymptotic value of the separation H is not the same 
as its value at the close encounter. We must also ask whether the eccentricity and inclination are 
affected in a similar way. In this situation e and i are not perfectly conserved adiabatic invariants 
(frequencies of the two superposed motions are not very strongly different from each other) but we 
probably will not make a huge mistake by assuming that they still are conserved on the approach 
and departure stages. Thus, we can draw the following picture of scattering for these orbits: as 
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Fig. 3. — Scattering probability P(Hq, H) as a function of the final orbital separation H for a fixed 
initial separation Hq. The probability is calculated by binning the outcomes of orbit integrations 
(histogram) and also analytically in Appendix A (solid line). Panels (a, b) display P(Ho,H) for 
fixed values of the initial eccentricity and inclination e = 2i = 5 and two values of the initial 
orbital separation, Hq = 2.05 and 4.85. Panels (c,d) display P(Hq,H) for fixed values of the initial 
dispersions of eccentricity and inclination, a e = &i = 5. Note the absence of scattered orbits near 
H = caused by the angular momentum exchange at large azimuthal distances (see text). 
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the planetesimal gets closer to the embryo the absolute value of H gets smaller as a result of the 
angular momentum exchange with the embryo, while the eccentricity and inclination are preserved 
(see Figure 4). At the point of closest approach, the planetesimal orbital parameters experience a 
quick variation as a result of close encounter with the embryo. After that, on the departure stage, 
the embryo slowly changes the angular momentum of the planetesimal such that \H\ increases. 

These variations of H at constant e and i do not contradict the conservation of the Jacobi 
constant (5) because the gravitational potential due to the embryo changes in the course of plan- 
etesimal approach and departure as well. In fact, using the conservation of Jacobi constant we 
derive a prescription for the relation between the semimajor axis difference H far from the embryo 
and its value at the moment of the close encounter, H'. To do this we note that at the moment 
of the closest approach the distance between the embryo and planetesimal in units of the embryo's 
Hill radius is of order \f e 2 + i 2 . Then the conservation of the Jacobi constant (5) tells us that 
H 2 — (H') 2 ~ (e 2 + i 2 ) -1 / 2 . This empirical expression may be inaccurate if one applies it for 
e, i < 1. This can be remedied using the following approximate prescription 

(H') 2 =H 2 - C ^ c~1.8, d~2, (27) 

V e 2 + i 2 + d 2 

where the numerical values of constants c and d were fixed by comparison with orbit integrations. 
Whenever the planetesimal disk has a distribution of eccentricities and inclinations we replace e 
and i in (27) with a e and cjj. The relationship (27) is certainly not very accurate. However, 
this transformation is satisfactory for our purposes since the complete separation of horseshoe and 
passing orbits we are assuming is a rather crude approximation anyway. 

To better understand the consequences and importance of this effect it is instructive to look at 
the distribution function of planetesimal scattering P(Ho,H) in the dispersion-dominated regime 
[by definition P(Hq, H)dH is the probability that a planetesimal initially at Hq is scattered into 
the interval (H, H + dH)\. We calculate this function analytically in Appendix A and represent it 
in the following form: 

P(H m- 1 {{Afl)2) + {Afl)AH AH-H H 

P{H(hH) ~2\nA (\AH\+d)* ' AH - H ~ H °> (28) 



where 

' - 

It \h\-\ 

d 



((a/*) 2 ; 

21nA 



(29) 



(Ah) and ((Ah) 2 } are the scattering coefficients (see Paper II), and In A is the usual Coulomb 
logarithm. To compute P(Hq, H) accounting for the distribution of planetesimal eccentricities and 
inclinations one simply needs to use the values of (Ah) and ((Ah) 2 ) averaged over this distribution 
(they are given in §4 of Paper II for the case of Gaussian distribution of e and i). 

One can also obtain this distribution function from numerical orbit integrations. To do this 
we perform Monte-Carlo simulations in two different ways: by fixing the absolute values of the 
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eccentricity and inclination of the planetesimals (only their epicyclic phases are chosen randomly), 
and by picking their orbital elements from the distribution (3) with fixed dispersions. In both cases 
the initial difference of semimajor axes Hq is the same for each simulation. In Figure 3a,b we plot 
P(Ho, H) for fixed e = 2i = 5 and in Figure 3c, d we do this for fixed a e = 2<Ji = 5. The analytical 
distribution given by (28) is also plotted for each case. One can see that it generally agrees very well 
with the numerical results. The large variations of P(Hq,H) in the outer wings of the simulated 
distributions are caused by statistical noise. When \H — Hq\ < 1 the agreement between theory 
and simulations is quite remarkable even for the distributions averaged over e,i (Figure 3c, d). 

However, one can also immediately notice one important feature of the simulated P(Hq,H): 
there are no orbits scattered in a region with a width ~ Rh around the embryo's location, and 
this is in contrast with the analytical result (28) which exhibits no such feature. A similar "gap" 
in the distribution of scattered orbits can also be noticed in the numerical calculations performed 
by Greenzweig & Lissauer (1990) and Ohtsuki & Tanaka (2002). 

From our previous discussion the reason for the appearance of this gap becomes quite apparent. 
It is the gravitational interaction of the particles on passing orbits located close to the horseshoe 
region with the embryo which drives them away from its orbit, corresponding to H = 0. If we were 
to take the same probability distributions not far away from the embryo in the azimuthal direction 
but immediately after the scattering (within ~ Rh from the embryo along the azimuthal direction) 
we would not find such a "gap" . However, as the differential shear slowly increases the azimuthal 
separation of the interacting bodies their mutual angular momentum exchange leads to a gradual 
increase of \H\ similar to that happening on the horseshoe orbits. As a result, a conspicuous "gap" 
appears in P(Hq, H) near H = 0. This process is illustrated in Figure 4. 

This effect has important consequences for the calculation of scattering on passing orbits with 
radial separations ~ Rh (or \H\ ~ 1). The difference between H far from the embryo (which we 
take as one of the planetesimal coordinates) and its value at the moment of the close encounter 
introduces changes into the computation of scattering coefficients. Indeed, if before and after the 
encounter the values of planetesimal semimajor axis are Hq and H, and at the point of closest 
approach it has the value of H', we can write that (see Figure 4) 

H'(H ) + Ah'[H'(H )] = H'(H + Ah), (30) 

where Ah' and Ah are changes of H' and H correspondingly. Our calculation of scattering coeffi- 
cients in §4 of Paper II gives us only the value of Ah' (change of the planetesimal semimajor axis at 
the closest approach). Here however we are interested in Ah = H — Hq. To relate them we assume 
that both changes are small which is always a good approximation in the dispersion-dominated 
regime. Then, expanding the r.h.s. of (30) up to the second order in Ah we find that 

A ^°> -(H)"' ^' - (if ) " 3 ^ + °«^> 3 >- 

(f)H'\ ~ 2 
^-J (Ah') 2 + 0((Ah'f). (31) 
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Fig. 4. — Schematic illustration of planetesimal scattering on passing orbits near the horseshoe- 
passing orbits boundary. The thick solid line in the left panel represents the track of the guiding 
center of some planetesimal initially located at Hq. As a result of angular momentum exchange 
with the embryo its semimajor axis shifts to H' prior to scattering; then in the course of a close 
encounter it is changed by Ah'. Subsequent interaction with the embryo on the departure stage 
increases the semimajor axis from H'+Ah' to H. Note that Ah = H — Ho is not equal to Ah'. Thin 
solid lines show the trajectories of planetesimals with the same Hq but different epicyclic phases. 
The panels on the right illustrate the origin of the gap in the probability distribution of scattering 
by plotting snapshots of the semimajor axes distribution of planetesimals P(Hq, H) initially (where 
it has (5-function form, at point A), immediately after the encounter (at point B), and long after 
(points C and C coincide in the closed box). 
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At the same time our assumption of approximate adiabatic conservation of e and i implies that the 
changes of these quantities evaluated far from the embryo are the same as they are in the course of 
an encounter near the embryo. It also means that in calculating scattering coefficients we should 
use the values of eccentricity and inclination characterizing the planetesimal orbit far from the 
embryo. 

The relationships represented by equations (31) are utilized when we calculate scattering co- 
efficients in equations (25)- (26). 



3.3.1. Boundary conditions for the scattering on passing orbits. 



The simplest way to derive the boundary conditions for equations (25)-(26) is to consider 
the Fokker-Planck equation for the evolution of the distribution function of planetesimals f(T) 
caused by embryo-planetesimal encounters. Here T is a set of orbital elements characterizing the 
planetesimal orbital state: Fi = H, (I^,^) = e (vector eccentricity) and (r^Ts) = i (vector 
inclination). Following conventional wisdom (Lifshitz Sz Pitaevskii 1981; Binney & Tremaine 1987) 
we can write the evolution equation of / caused by weak encounters: 



df 
dt 



= -2\A\ 



d 2 



1=1 1=1 J 



(flHKATiATj)) 



dF 



where the flux in the i-th. direction is given by 



2\A\ 



f\H\{ATi 



TUT* 02) 



(33) 



Averaging of quantities like (Ar^) is performed here only over the possible outcomes of scattering 
and not over dedi. The factor 2|^4||i/| takes care of the linear shear velocity in the planetesimal 
disk (changes of various quantities caused by the scattering are assumed to be per encounter and 
not per unit of time). Note that we could have obtained equations (25)-(26) directly from (32). 

Particles on passing orbits do not penetrate into the region of horseshoe orbits (as a result of 
our presumed complete separation of these two types of orbits). Thus the component of the flux 
Fh must vanish at the boundaries H = hh s and H = —hh s : 



[f\H\(AHA(e 2 ))] 



1 d 



2d(i 2 ) L 



f\H\(AHA(i 2 )) 



l^[m\((AHf)} 



at H = ±h hs . 



(34) 



Our next step is to multiply condition (34) by d(e 2 )d(i 2 ), e 2 d(e 2 )d(i 2 ) and i 2 d(e 2 )d(i 2 ) and 
integrate it over e 2 -i 2 space. We obtain as a result (integrating by parts where needed) that 



N\H\(Ah)-^[N\H\((Ah) 2 ) 



0. 



(35) 
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N\H\{e 2 Ah) - \ A [N\H\(e 2 (Ahf)] + liY|^|(A(e 2 )A^) = 0, 



(36) 



N\H\{i 2 Ah) 



1 d 
2~dH 



1 



N\H\{i 2 (A~h) 2 )\ + -JV|ir|<A(?)A/i> =0 at H = ±h hs , (37) 



where now (...) means integration not only over the probability distribution of scattering but also 
over the e 2 -i 2 space (and we use the more familiar notation {Ah) instead of {AH), etc.). In getting 
these results we have taken into account that all scattering coefficients vanish at the boundaries of 
the velocity space (which are at infinity). 

As we discussed in Paper II terms like ((Ae) 2 A/i) are third order in perturbations and can be 
neglected. Thus we can approximate {A(e 2 )Ah) with 2((e • Ae)Ah). As a result we find that 



N\H\{A~h) 



1 d 
2dH 



N\H\{(A~h) 2 ) 
1 d 



= 0, 
iY|#|(e 2 (A/i) 2 ) 



N\H\{(e 2 +e-Ae)A~h) ^ 
N\H\{(1 2 + i • Ai)Ah) - ± A \N\H\(J 2 (Ah) 2 ) 



0, 



at H = ±h 



hs- 



(38) 
(39) 
(40) 



2dH V 

Note the specific combination e 2 + e • Ae, not e 2 + 2e • Ae as in equation (26) — this is important 
for the conservation of the Jacobi constant. 

Do equations (25)-(26) conserve the integrated Jacobi constant (6) and the total number 
of planetesimals (7)? If we integrate (25) over dH (which would give us the total number of 
planetesimals in the region of interest), we find that condition (38) ensures the conservation of 
N tot . The same is true for the integrated Jacobi constant: when we substitute equations (25)- 
(26) into the definition (6) we find (after some cumbersome but straightforward calculations) that 
conditions (38)- (40) guarantee the conservation of J tot up to the second order in the perturbed 
quantities. Thus, expressions (38)-(40) provide a desired set of self-consistent boundary conditions 
for the equations (25)- (26). 



3.3.2. Comparison with numerical orbit integrations. 

To check our predictions about the behavior of the planetesimal disk properties derived in 
the previous sections we have performed a set of numerical simulations. We have integrated the 
orbits of test particles, starting at large azimuthal separation from the embryo, and observe the 
changes of their orbital parameters as they experience gravitational interactions with the planetary 
embryo. Initial orbital parameters are chosen randomly from the distribution of eccentricities 
and inclinations (3). The semimajor axes of particles are assumed to be uniformly distributed 
within some radial interval around the embryo. The width of this interval is large enough that 
boundary effects are not important. We typically calculate about 2 x 10 5 different orbits to reduce 
statistical noise. Each orbit experiences several hundred passages past the embryo during which 
its orbital elements change. Between the conjunctions we randomize the epicyclic phases of the 
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planetesimals to mimic the effect of planetesimal-planetesimal interaction which is assumed to 
destroy any resonances in the system (the hypothesis of molecular chaos). However the absolute 
values of eccentricity and inclination are not affected by this procedure. In the course of the 
integration embryo is assumed to be immobile. 

The number of consecutive passages is dictated by the condition that the system evolves for at 
least one dynamical time within region of phase space for which close encounters can occur \H\ < a e . 
In the shear-dominated regime this corresponds to a single passage at radial separation ~ Rh (i.e. 
dynamical time td yn ~ tsyn — synodic period at H = 1) because scattering is strong in this case. 
At the same time equations (25)- (26) and analytical expressions for the scattering coefficients in 
Paper II show that in the dispersion-dominated regime the dynamical time td yn ~ t syn a^ when 
a e ~ (jj. This timescale can easily take several hundred orbital passages which require that one 
follow planetesimal orbital evolution for rather a long time. 

In Appendix B we describe our treatment of the intermediate (between shear- and dispersion- 
dominated) velocity regime and make comparisons with the results of orbit integrations. Here 
we perform this procedure for the dispersion-dominated regime. In Figures 5 and 6 we display 
time sequences of profiles of horizontal and vertical velocity dispersions [panels (a,b) and (c,d) 
correspondingly], and dimensionless surface density normalized by its value at infinity [panels (e,j)]. 
We show both the results of orbit integrations (left row of panels) and our analytical predictions 
(right row of panels) . Disk evolution in the region of the horseshoe orbits is described using equation 
(24). In the region of passing orbits [outside of the interval (—hh s -,hh s ), see §3.1] disk evolution is 
governed by partial differential evolution equations (25)-(26). We solve them using fully implicit 
scheme (Press et al. 1988) with the boundary conditions (38)-(40) imposed at H = —h^g and 
H = hh s and scattering coefficients computed analytically in Paper II. The conversion (27) is used 
throughout the calculation and the boundary of the horseshoe region is described by formula (22). 
For the factor A entering the Coulomb logarithm in the expressions for the scattering coefficients 
we use the following prescription 4 : 

A = <7i(2<7g + erf). (41) 

Constant coefficients in this formula are roughly fixed using comparison with orbit integrations 
(but our final results depend on them very weakly). 

Figure 5 shows the disk evolution in the case when the initial planetesimal dispersions of 
eccentricity and inclination are both equal to 3. One can see an excellent qualitative agreement 
between the results of orbit integrations and analytical theory. All the features of the spatial 
distributions of disk quantities are well reproduced by the solutions of the analytical equations. 
There is also a reasonable degree of quantitative agreement between them although there are also 
some minor differences: the analytical equations predict somewhat faster evolution of &i and a 



4 To avoid problems at a e ,ai < 1 we use instead of InA the more accurate expression (1/2) ln(l + A 2 ) (see Binney 
& Tremaine 1987). 
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Fig. 5. — Evolution of the planetesimal disk properties driven by the presence of a massive pro- 
toplanetary embryo. Initially a e = cr, = 3. The plots contain numerical (left row) and analytical 
(right row) time sequences of profiles of a e (a,b), di (c,d), and dimensionless surface density normal- 
ized by its value at infinity (e,j). Curves of different colors represent profiles measured at specific 
moments of time normalized by the synodic period at a separation of H = 1 (which are shown 
in panels (e) and (f) by corresponding color coding; they are slightly different for numerical and 
analytical curves). 
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Fig. 6. — The same as Figure 5 but for initial a e = 4, <7j = 2. 
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slightly smaller radial extent of the excited region than orbit integrations do (this might be caused 
by the effect of the distant encounters which theory can not reproduce very well, see Appendix 
B). In Figure 6 we display the same comparison of numerical and analytical results but now for 
the case of a e = 4, di = 2. This corresponds to the ratio di/a e = 0.5 which is within the range 
0.45 — 0.55 thought to be realized in homogeneous Keplerian planetesimal disks (Ida et al. 1993; 
Stewart & Ida 2000). One can see again a very good consistency between the two approaches, 
especially for the evolution of N and <jj. There is a slight discrepancy between the predictions of 
two methods for the evolution of a e near the horseshoe-passing orbits boundary. However, since 
equations (25)-(26) preserve the Jacobi constant one can be sure that the evolution of the surface 
density (which is driven by the redistribution of the angular momentum in the planetesimal disk) 
is always consistent with the evolution of epicyclic energy independent of the details of the spatial 
distributions of various quantities. 

Planetesimals with the semimajor axes close to the embryo's orbit perform horseshoe motion 
(see §3). As a result, surface density of guiding centers in this region stays constant in time as can 
be seen in Figures 5 and 6. Note that our orbit integrations do not exhibit the concentration of 
planetesimals at the embryo's semimajor axis which is present in numerical calculations of Tanaka 
& Ida (1997). In their case this effect might have been caused by the analytical simplifications 
which were employed in Tanaka & Ida (1996, 1997) to speed up orbit integrations. 

Figures 5 and 6 show that the embryo does not clear a complete gap on a dynamical timescale 
in the dispersion-dominated regime, although it does produce a significant depression in the surface 
density. Note that we plot the surface density of guiding centers of planetesimals N in Figures 5 
and 6 rather than the instanteneous surface density. The depression formed in the instanteneous 
surface density at the location of the embryo is weaker than the gap in N because epicyclic motion 
allows planetesimals to penetrate inside the depression of N. The absence of a clear gap is due to 
the weakness of the individual embryo-planetesimal scattering in the dispersion-dominated regime. 
In the shear-dominated disk, scattering is much stronger and a gap is cleared by the embryo on a 
rather short timescale (see Paper I and Appendix B). It is also obvious that the evolution of the disk 
surface density is accompanied by a considerable change in kinematic properties of planetesimal 
population, as required by the conservation of Jacobi constant. Thus, the dynamical evolution of 
the disk is an important ingredient of the gap formation process which justifies the need for the 
self-consistent theory such as the one presented here. 

The agreement between the analytical and numerical results in the dispersion-dominated 
regime is in general pretty good given the approximate nature of the theory and the small number 
of fitting parameters we are using — essentially only the constants in (20), (27), and (41) are free 
parameters, and it turns out that the solutions of equations (25)-(26) depend on their particular 
choice only very weakly. We believe that minor discrepancies between the outcomes of analytical 
and numerical approaches can be accounted for by going to the next order in 1/lnA and properly 
including the effects of the distant encounters. We expect our analytical formulation to be even 
more accurate for larger values of a e , &i but we have not investigated these because the compu- 
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tational requirements of the numerical simulations become prohibitive. We postpone the detailed 
exploration of these subjects for the future. 



4. Embryo's accretion rate. 

The accretion rate of planetesimals by the embryo M is another observable (in addition to the 
spatial distributions of N, a e , and cjj) which can be computed both analytically and numerically 
and used as a check of our calculations. In addition, it is a very important quantity by itself 
for any realistic modelling of planet formation. A lot of work has been devoted to the study of 
the planetary accretion in homogeneous planetesimal disks in the context of the planetary mass 
growth (Greenzweig & Lissauer 1990, 1992) and the origin of planetary spins and obliquities (Dones 
& Tremaine 1993). The results of these authors are not directly applicable to the inhomogeneous 
disks that are our primary focus but can be used as limiting cases to check our analytical predictions 
in the more general case of nonuniform distribution of planetesimals. 

In Appendix C we calculate M in different velocity regimes. In the dispersion-dominated case 
such a computation is made possible by the use of the two-body approximation and we find that 



M 



m- 



, . \H\dH 



(i/')7(25 e 2 ) 



(42) 
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H'/dj 
2^2 



(43) 



R e is the embryo's physical radius, a e is its semimajor axis, p = R e /Rn, Iq is a modified Bessel 
function of order zero, and H'{H) is given by (27). Integration over dH in (42) should exclude the 
region (—hh s ,hh s ) corresponding to the horseshoe orbits. 

In the shear-dominated regime we use simple scaling arguments to fix the accretion rate de- 
pendence on the physical variables of the system and find that 



M « 5N mst (H coU ) 



mQR e R H 



{di)a 



(44) 



where N mst (H co u) is the surface density of planetesimals on orbits leading to collisions with the 
embryo (at H = H co u rs 1.4, see Appendix C), defined by equation (CIO), and (<7j) is a measure of 
cjj in the inhomogeneous disk obtained by averaging Gi over the interval (—2Rh, 2Rh) [from where 
most of the planetesimals get accreted in the shear-dominated regime, see Petit & Henon (1986)]. 
In Appendix B we demonstrate how to use these results to cover the regime of intermediate velocity 
dispersion. 
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Fig. 7. — Dimensionless accretion rate m = M/(f2i?g£ ) (-^e is the radius of the planet, E is the 
surface mass density of the disk far from the embryo) as a function of time (expressed in synodic 
periods t syn at the separation of 1 Rh) for different values of a e and &i. Results of orbit integrations 
are shown by thick curve, and the dotted line shows analytical predictions. Accretion rates for the 
dispersion-dominated regime are presented in panels (a) and (6), while panels (c) and (d) display 
rh for intermediate velocity regime. 

Figure 7 shows the comparison of our analytical predictions with the results of orbit integra- 
tions. We plot the dimensionless accretion rate rh = M/(QR1T,q) (i? e is the radius of the planet, 
So = mN(oo)/a^, is the dimensional surface mass density of the disk far from the embryo) as a 
function of time normalized by the synodic period t syn corresponding to a separation of 1 Rh be- 
tween the orbits of the embryo and planetesimal. Numerical accretion rates were computed during 
the calculation of the spatial and temporal evolution of various disk quantities described in §3.3.2. 
They are displayed by the thick solid line in Figure 7. Analytical results [using an interpolation 
given by formulae (B5) and (B6) where needed] are shown by a dotted line. One can see that the 
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agreement is reasonably good even in the intermediate velocity regime 5 (Figure 7c, d). The theo- 
retical results in the dispersion-dominated regime (Figure 7a,b) do not depend on the interpolation 
and their good agreement with the numerical accretion rates confirms the validity of the analytical 
approach developed in Appendix C. We will be using this prescription for the embryo's accretion 
rate in upcoming work when dealing with the evolution of the planetesimal disk coupled to the 
embryo's growth. 

5. Direction of the embryo-planetesimal interaction. 

From the results of the orbit integrations described in §3.3.2 one can see that the embryo- 
planetesimal interaction leads to a decrease in the surface density of planetesimals near the embryo 
in the region of passing orbits. The width of the surface density depression is typically of the order 
of the planetesimal eccentricity dispersion a e . At the same time the distribution of planetesimals 
on horseshoe orbits stays unaffected by the embryo. 

This is very similar to the behavior of the surface density in the shear-dominated regime, 
as studied in Paper I. In that case the reason for such behavior was very clear: scattering of a 
planetesimal initially on a circular orbit can only increase its eccentricity and inclination which 
leads to the increase of \H\ (a result of Jacobi constant conservation) and, consequently, to the 
repulsion of planetesimal orbits. In Paper I planetesimals were kept on circular orbits at all times, 
and a gap with a width of several Rh was carved near the embryo's orbit. 

In the kinematically hot regime the reasoning is not so simple. From equation (28) and Figure 
3 one can see that scattered planetesimals can have A\H\ < as well as A|if| > for a fixed H 
depending on their eccentricities, inclinations, and epicyclic phases. Moreover, from our calculation 
of the scattering coefficients in the dispersion-dominated regime (Paper II, see also the analogous 
result in Ida et al. 2000) we know that (Ah) < 0. Using this result Ida et al. (2000) suggested that 
in the kinematically hot planetesimal disk embryo should attract planetesimal orbits rather than 
repel them. The outcome of such an interaction would be a growth of planetesimal surface density 
at the embryo's location. 

However, this conclusions is misleading. What really determines the behavior of the planetes- 
imal surface density near the embryo and the nature of the embryo-planetesimal interaction is the 
direction of the flux of planetesimals. The embryo repels planetesimals and carves a gap when 
this flux is directed away from the embryo; it attracts planetesimals and increases their nearby 
surface density when this flux is directed towards the embryo. The condition (Ah) < by itself 
cannot guarantee that the planetesimal flux is directed towards the embryo's orbit because these 



5 Note that the specific form of interpolation function (B6) was chosen to fit several sets of a e and a e , not only 
those displayed in Figure 7. In all studied cases numerical and analytical results agree within (10 — 30)%. 
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Fig. 8. — Contour plot of the dimensionless planetesimal flux ttT /(\A\^ e ) as a function of ai/a e 
- ratio of the vertical to horizontal velocity dispersions, and H/a e — radial distance from the 
embryo scaled by the RMS epicyclic radius. Note that this flux is always positive for positive H 
meaning that the embryo repels planetesimals. The flux is an odd function of H. 

quantities are not simply related 6 . Another complication precluding the use of the sign of (Ah) as 
an indicator of planetesimal disk evolution is that AH < (for H > 0) does not always mean that 
planetesimal orbit gets closer to the embryo. If for example AH < —2H then the absolute value of 
the semimajor axis separation after the encounter is larger than it was before, which is indicative 
of repulsion rather than attraction. 

Using the method described in Appendix A of Paper I, we can calculate the outward flux 
of planetesimals T driven by the embryo-planetesimal interactions. We evaluate this flux in a 
homogeneous disk. For the general scattering probability function P(Hq, AH) [P(Hq, AH)dAH is 
the probability for a planetesimal initially at Hq to have AH in the range (AH, AH + dAH)] this 



6 One can imagine a scattering probability function which transports most of the planetesimals away from the 
embryo shifting them by a small AH > 0; at the same time a small fraction of planetesimals can be scattered towards 
the embryo and have large negative values of AH. With properly chosen weights one can always ensure that {Ah} < 
while the bulk of material moves away from the embryo forming a gap around its orbit. 
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flux can be written as (see Papers I and II) 
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Differentiating this expression w.r.t. H gives us 
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Since the planetesimal flux in a homogeneous disk has to vanish far from the embryo we can 
integrate (46) to obtain 
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(47) 



This expression [or (45)] is valid for an arbitrary function P(Hq, AH). In principle we can use 
this general form but it is more convenient to apply one of the most important properties of the 
interaction in the dispersion-dominated regime: that scattering is dominated by weak encounters, 
meaning that we can expand in (47) the integral in brackets in a Taylor series in AH. As a result, 
after some intermediate transformations we find that planetesimal flux in the initially homogeneous 
disk is given by 
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(48) 



This expression could have also been obtained directly from equation (25) or (33). It is clear from 
the equation (48) why the sign of (Ah) cannot tell us the direction of the planetesimal evolution: 
the presence of additional term with the derivative of \H\((Ah) 2 ) in r.h.s. of (48) which cannot in 
general be related 7 to (Ah) means that there is no direct relation between the signs of F(H) and 
(Ah). 

Using the analytical formulae for (Ah) and ((Ah) 2 ) derived in Paper II, we can calculate this 
flux in the Fokker-Planck approximation as a function of distance from the embryo and the ratio 
of the vertical to horizontal velocity dispersions. The result of such a calculation is presented in 
Figure 8. The shape of the contours is universal in the coordinates H/a e ,di/a e , i.e. it does not 
depend on the embryo mass, planetesimal mass, or other parameters. In plotting T we have also 



7 A relation between the first and second order diffusion coefficients exists if a homogeneous distribution is a 
steady-state solution of the Fokker-Planck equation see Lifshitz & Pitaevskii (1981). In our case this is not true — 
embryo's gravity tends to make disk inhomogeneous by carving out a gap. 
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neglected the effect of interaction at large azimuthal distances near the horseshoe region boundary 
described in §3.3 (it can be easily taken into account and will only slightly shift the contours near 
H ~ 1.) One can easily see from Figure 8 that T > when H > for all values of H/a e and (Ji/a e 
which we explored (and T < for H < by symmetry). This implies that the embryo tends to 
repel planetesimals by driving a flux directed away from its orbit. This conclusion disagrees with 
the predictions based on the sign of (Ah) (Ida et al. 2000). The planetesimal flux vanishes as H/a e 
grows but distant encounters (not included in Figure 8) will dominate at large H/a e and they also 
always drive planetesimals away from the embryo (Goldreich Sz Tremaine 1982; Petit & Henon 
1987a). Thus, it appears that repulsion is the only possible outcome of the embryo-planetesimal 
gravitational interaction even in the dispersion-dominated velocity regime. 

Figure 8 seems to imply that the planetesimal surface density will always grow because T 
decreases with H. This is not the case because the boundary condition (38) forces T = at the 
horseshoe-passing orbit boundary H = hh s . This means that there is actually a transition zone 
near H = ±hh s in which flux steeply increases from to its value given by (48). The width of 
this zone is infinitesimal initially but it rapidly expands as the planetesimal disk evolves under the 
action of the embryo's gravity. Within the transition region T is positive and increases with H 
meaning that the planetesimal surface density decays there. This is exactly what was observed in 
the analytical calculations and numerical orbit integrations in §3.3.2. The slowdown of the disk 
evolution obvious from Figures 5 and 6 is caused not only by the growth of epicyclic energy of 
the disk (which lengthens the disk dynamical timescale) but also by the gradual expansion of the 
transition zone leading to a decrease in the planetesimal flux gradient and slower disk evolution. 



6. Discussion and conclusions. 

We have studied the embryo-planetesimal interaction in the gravitational field of a central star. 
Two different cases were explored: when the interaction between the embryo and planetesimals 
occurs in the shear-dominated and in the dispersion-dominated regimes. The treatment of the first 
case parallels that of Paper I but is complicated by the fact that now we explicitly include the 
evolution of not only the surface density but also the eccentricity and inclination of planetesimals 
in the disk. 

Our study of the dispersion-dominated regime relies on the methods of kinetic theory, and it 
uses many of the results obtained in Paper II. However our present treatment is more refined since 
the description of embryo-planetesimal scattering requires clarifying many details which were not 
important for the planetesimal-planetesimal interactions. In particular we have to study not only 
passing but also horseshoe orbits of planetesimals to determine the spatial distribution of the disk 
properties. To do this we propose a condition which separates the horseshoe and passing orbits and 
check its viability using numerical orbit integrations. Angular momentum exchange between the 
embryo and planetesimal long before and after their closest approach turns out to be important for 
the scattering on passing orbits near the horseshoe boundary. We illustrate this point by comparing 
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the analytical scattering probability function with the one obtained from numerical integrations. 
A simple method to account for this effect in our Fokker-Planck approach is proposed. Taken 
together all these refinements are shown to provide rather good agreement with the results of the 
numerical orbit integrations. Thus we hope to have grasped the most important features of the 
embryo-planetesimal interaction by our theoretical approach. 

This does not mean that our treatment of the embryo-planetesimal interaction is complete. 
We have only focussed on the most important, dominant effects, and there is certainly room for 
additional refinements, which would farther improve the agreement with numerical results. The 
calculation of the scattering coefficients to the next order in 1/ln A would provide us with subdom- 
inant corrections which might be important for modest values of a e and o"j. One can certainly do 
a better job in treating distant encounters, calculating various coefficients entering formulae (20), 
(27), etc. or the interpolating functions of Appendix B, using a larger set of numerical orbit inte- 
grations. On a somewhat deeper level, one can try to come up with a more sophisticated treatment 
of the horseshoe-passing orbits separation (instead of the complete spatial separation of these two 
types of orbits assumed in this paper) . Our purely deterministic treatment of the shear-dominated 
regime can also be improved, which would ameliorate the comparison with numerical results in 
the intermediate velocity regime (see Appendix B). The recoil of the embryo and the excitation of 
the embryo's eccentricity and inclination in the course of the planetesimal gravitational scattering 
can be important in some applications, such as the embryo's migration (Tanaka & Ida 1999) or 
its interaction with the embryos nearby (Tanaka & Ida 1997). Our treatment relies on the use of 
the Schwarzschild velocity distribution which was demonstrated to be applicable in the dispersion- 
dominated regime (Greenzweig & Lissauer 1992), but the deviations from this assumption could be 
important e.g. in the intermediate velocity regime, and this subject can also be pursued farther. 
All these refinements would better the quantitative agreement between the analytical theory and 
numerical results. But reasonably good accord is provided even by our basic treatment developed 
here, especially in the dispersion-dominated case where the assumptions we make are the most 
justifiable. 

We also dwell upon the question of the direction of the embryo-planetesimal interaction, namely 
whether it leads to the repulsion of planetesimal orbits from the embryo or to their attraction. The 
latter outcome has been favored in some scenarios (Ida et al. 2000) and is based on the fact that in 
the dispersion-dominated regime embryo on average tends to attract planetesimal semimajor axes 
toward its orbit. We demonstrate, however, that the average change of planetesimal semimajor 
axes cannot serve a standard for determining the direction of the embryo-planetesimal interaction 
because the transport of the angular momentum (associated with the changes in semimajor axes) 
is not the same as the bulk motion of the disk material. We propose our own criterion for judging 
the embryo-planetesimal scattering outcome, and show that the embryo always repels planetesimals 
rather than attracts them in an initially homogeneous disk thereby carving out a gap in the dis- 
tribution of the planetesimal semimajor axes [which is in contrast to claims by Ida et al. (2000)]. 
Our own numerical results support this conclusion. 
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The combination of our results for planetesimal-planetesimal gravitational scattering presented 
in Paper II and the theory for the embryo-planetesimal interaction developed here now allows us 
to study the planetesimal disk dynamics with these processes operating simultaneously. It also 
provides a framework to which various refinements and additional physical mechanisms can be 
naturally added. Our results are not restricted in applications to the problem of the formation of 
planetary systems but can also be used for studying their present day evolution, e.g. the dynamics 
of asteroid and Kuiper belts affected by massive bodies inside or near them. Our results for the 
accretion rate of massive body immersed in inhomogeneous planetesimal disk (§4 Sz Appendix C) 
allow us also to include self-consistently the embryo's mass growth into consideration. We will 
examine protoplanetary disk evolution with all these effects working together in a future work 
(Rafikov 2002b). 

I am grateful to my advisor, Scott Tremaine, for his patient guidance and valuable advice. 
The financial support of this work by the Charlotte Elizabeth Procter Fellowship and NASA grant 
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A. Scattering probability function in the dispersion-dominated regime. 

In this appendix we calculate the scattering distribution function P(H, AH) in the dispersion- 
dominated regime. This function is defined such that P(H, AH)d(AH) is the probability that two 
bodies initially separated by H will have a change of H in the interval (AH, AH + d(AH)) during 
an encounter. This probability will initially be considered as a function of also e and i. Then we 
will average it over the distribution of e and i and obtain the averaged value of P(H, AH) as a 
function of a e and 0%. We will extensively use the results and notation of §4 of Paper II which 
should be consulted for further details and clarifications. 

To do this we will employ the two-body approximation. We first relate AH to the change of 
the relative velocity in the y-direction Av y using formula AH = 2Av v /(QRh)- Then we switch to 
the coordinates I (impact parameter) and <p (angle in the plane perpendicular to the direction of 
approach velocity 8 ) to characterize the trajectory of the bodies. Using these coordinates one can 
express Av y as (Binney & Tremaine 1987) 

Av v = -At;,|^ + At; ± ^i^cos0, (Al) 

where v x ,v y ,v z and v are different components of the planetesimal approach velocity and its mag- 
nitude given by equations (58) and (60) of Paper II; expressions for Av^ and Av± are given by (67) 
ibid. This allows us to express AH as 

AH = ^T77\2 l-H + 2 (l/h) n cos 4>) • (A2) 
1 + (//to) 

where v\ = e 2 + i 2 — H 2 (vf is always positive because \H\ < e for the kind of scattering we 
consider here) and Iq = G(mi + m2)/v 2 = Rh(^th /v) 2 (when minimum approach distance between 
interacting bodies < Iq large-angle two-body scattering takes place). 

For a fixed AH this equation can be rewritten as an equation for u = (I /lo) cos <fi and w = 
(/// o )sin0: 



/ 25i \ 2 2 _ 2 2 _ (AH - AH.)(AH + - AH) 
{ U -Ah) +W - R > R ~ AJP ' 



(A3) 



where 



AH_ = -H - ^H 2 +4vf, AH + = -H + yj ' H 2 + \v\. (A4) 
Clearly this is only possible if 

AH- <AH< AH+. (A5) 



The origin from which this angle is counted is not important here. 
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Thus, the difference in planetesimal semimajor axes cannot increase by more than AH+ or decrease 
by more than AH- during the scattering. Equation (A3) is simply the equation of a shifted circle 
with radius R in coordinates u, w. 

To calculate the probability distribution function P(H, AH) we use the following identity: 



£P(H,AH'<z) 



AH < 0, 



P(H, AH) = { , . . Z=AH (A6) 



-A 

dz 



P(H, AH' > z) 



AH > 0, 

AH' 



where P(H, AH' < z) and P(H, AH' > z) are complementary cumulative probability distributions. 
We calculate them in the following way. 

It is clear that the probability P(H, S) of some outcome of scattering is equal to the ratio of 
the area of the region of the space of epicyclic phases t — to corresponding to the outcome manifold 
S and 4-7T 2 . Using the conversion between the integration over r — to and / — <p represented by 
equation (61) of Paper II we can write this as 

P(H, S)=[*^ = -\~ 1= I IdA (A7) 

K ' > J 4vr 2 3tt R% QR h i\H\Ve 2 - H 2 J 2vr V ; 

s ' s 

Thus we can always transform our problem into the calculation of the area of region S in I — (j> 
coordinates. Here we are interested in the area of the region given by AH < z for z < and by 
AH > z for z > 0. 

One can easily see that both conditions correspond to the inner part of a circle given by 
equation (A3). Thus the area of the region in I — <f> phase space bounded by these conditions is 
always represented by a single formula (no matter whether z is positive or negative): 



/ 



^ = ^ = g( i -Ag z KAg ± - £ ) 

2vr 2vr 2 z 2 v ' 

s 



Substituting this result into (A7) and then into (A6) we obtain (using the definitions of v and Iq) 
that 

4 4v 2 - HAH 1 

P(H, AH) = i a 7 = tt^ (A9) 

3vr \AH\* i\HWe 2 -H 2 [e 2 + i 2 - {3/4)H 2 ] 3/ 

if AH- < AH < AH + , and P(H, AH) = otherwise. Recalling the definition of v\ and using the 
definitions of scattering coefficients (Ah) TjUI and ((Ah) 2 ) T ^ given by equations (73), (74), and (77) 
of Paper II we can rewrite (A9) in a more appealing form: 

P(H,e,i,AH) - — ^ , (A10) 

where the dependence on e and i is hidden in (Ah) T>(0 and ((Ah) 2 ) TjUJ . 
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The expression for the differential probability given by equation (A10) is most accurate for large 
angle scattering (i.e. for AH > 1) because this means very small impact parameters implying that 
the two-body scattering assumption is very good. This allows one to use this probability function to 
study strong scattering in the Hill approximation. On the other hand when AH — > the expression 
(A10) clearly diverges. This divergence is unphysical because the two-body assumption becomes 
very bad there. Small AH means weak scattering which only occurs for large impact parameters 
I. In real planetesimal disks I is always restricted to be < iRjj- Thus our formula (A10) is only 
applicable for \AH\ > AH m i n , where we find AH m i n by setting I zRh 3> lo in equation (A2): 

Ai^„^(^y«i. (ah) 



L min 

I 



\ V J 



When \AH\ < AH m i n the scaling provided by (A10) is not applicable and the probability 
distribution of scattering assumes some different form. However, in most cases we need not worry 
about the effects of this indeterminacy for small AH because usually we need only the integrated 
characteristics of the probability function [e.g. in calculation of some scattering coefficients in Paper 
II which could be done using P(H, AH), see below]. Integrals of P(H, AH) are subject to strong 
cancellation effects near AH = and this naturally removes the problem. Similar situation was 
discussed by Goodman (1983, 1985) in the context of stellar gravitational encounters in globular 
clusters. 

Using (A9) we can provide an alternative derivation of (Ah) TjUl . Multiplying r.h.s. of (A9) by 
AH and integrating from AH- to —AH m i n and from AH m i n to AH + one obtains that 9 

, ~, 4 h\n(\AHAAH + /AH 2 in ) 

(Ah)^ = - — U 1 +/ - mm> (A12) 

37r i\H\y/e 2 -H 2 [e 2 + i 2 - (3/4)H 2 ] S/ 

Using (A4) and (All) we find that 



which means that expression (A12) coincides with equation (73) of Paper II. 

For some purposes it will still be necessary to use the distribution probability (A10) itself rather 
than its integrals. For these occasions we can remedy the divergence in (A10) by introducing some 
cutoff distance d ~ AH m i n into the denominator of (A10). One can do this in a variety of ways, 
one of the simplest ones would be to take 

P{H M AH) = — ( , Ag | +d)3 • ( A14 ) 



9 Remember that we are using (Ah) and {(Ah) 2 ) instead of (AH) and {(AH) 2 ) for consistency with the notation 
of Paper II. 
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One can easily check that the normalization of the total probability to 1 would require that 

1/2 



d 



2' 



21nA 



(A15) 



We use this nondivergent form of (A14) in §3.3. 



To take into account the distribution of planetesimal eccentricities and inclinations we need 
to average the differential probability using Schwarzschild velocity distribution. Expression (A10) 
provides a natural way of doing this and we find as a result that 

p (H!i * AH)- 1 {{A ~ h)2) ± {Ajl)AH (AW) 
P(H,a e ,a t ,AH) - — ^ , (A16) 

where explicit analytic expressions for {Ah) and {(Ah) 2 ) are given by equations (86) and (87) 
of Paper II. In principle we should be bearing in mind the restriction (A5) when doing this last 
averaging. In practice this can only slightly affect the factor inside the logarithm if we restrict the 
range of AH by approximate condition 

AH- < AH < AH + , AH± = -H ± ^ H 2 + 4{a 2 + a 2 - H 2 ). (A17) 

Outside of this range simple form (A16) will give rather poor approximation. The nondivergent 
form of (A16) can be obtained in the same fashion as (A14) and (A15). 



B. Intermediate velocity regime and distant encounters. 

Equations (13)-(15) and (25)-(26) describe the evolution of the disk properties in the shear- 
and dispersion-dominated regimes correspondingly. Unfortunately the assumptions we have made 
while deriving them preclude us from using these equations directly in the intermediate regime, 
when 

of + cr? ~ 1. (Bl) 

Still we can try to provide an approximate description of evolution in this velocity regime by 
smoothly interpolating between the two extremes we believe we can describe well. Of course, one 
should not expect very good agreement of the theory interpolated into the intermediate velocity 
regime with the outcomes of orbit integrations. Qualitative concord would be enough for our 
purposes. 

Far from the embryo where the scattering due to close encounters is exponentially small (be- 
cause of the scarcity of planetesimals with large eccentricities) distant encounters dominate the 
evolution of the planetesimal disk because their effect decays only as a power law in the embryo- 
planetesimal separation. We attempt to take them into account by using the same interpolation 
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technique: we assume that when H/a e <C 1 close encounters dominate and scattering is described 
by formulae of §3.3. When H/a e » 1 distant encounters dominate and we use formulae of §3.2 
to describe their effect on the disk properties because distant scattering is similar to the shear- 
dominated scattering. 

To provide a smooth matching between the different velocity regimes and spatial scattering 
zones we use interpolating functions 0(x) and ^(x) with the following properties: 0(x),\E'(x) — ► 
as x — > and G(x),^(x) -> 1 as s > 1. Then we can describe evolution of some quantity F 
(N, &1,af) in terms of the corresponding shear- and dispersion-dominated extremes as 



9F ^,-2 ~2^ F 



+ [l-Q(a 2 e + a?)(l-*(H/a e ))] — 

disp 



(B2) 

shear 



The properties of G(x) and *(x) ensure that this formula reproduces correct limiting behaviors of 
disk evolution. 

We have found that the results of orbit integrations performed in the intermediate velocity 
regime can be satisfactorily described using the following form of interpolating function 0(x): 

40 



O(x) = exp 



(B3) 



X 1 /2( 5 + x )3_ 

For \I/(x) we use the following ad hoc prescription: 

*(x) = exp [-(8/x) 4 ] . (B4) 

Comparison of the disk evolution in the intermediate regime computed using equations (B2)- 
(B4) with that obtained from direct orbit integrations is shown in Figure 9. We display the results 
for a disk with initial a e = Gi = 1. One can see that the agreement between two methods is 
only qualitative indeed. The most probable reason for this is the use of the simple prescription 
(13)-(15) with a discrete probability function (12) to treat the shear-dominated regime. This also 
rules out the possibility of accurate treatment of the effects of distant encounters. In fact, as 
long as the random motion in the disk in not exactly zero scattering always has some dispersion 
of outcomes and this gives rise to the diffusive transport of planetesimals in the disk. Inclusion 
of this effect would widen and lower the surface density profile as well as remove spiky features 
caused by the discrete scattering function (12). We do not pursue this subject here and leave its 
detailed exploration for the future. Still, we have verified using several different choices of a e and 
<7j in the intermediate velocity regime that qualitative agreement is adequate: general features and 
trends appearing in the numerical results are always reproduced reasonably well by our analytical 
prescription. 

In Appendix C we describe the calculation of the planetesimal accretion by the embryo in 
the shear- and dispersion-dominated regimes. We again use the interpolation approach for the 
description of the accretion rate in the intermediate velocity regime: 

M = Ii{d 2 e + af)M disp + [1 - n(a e 2 + af)] M shear , (B5) 
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Fig. 9. — The same as Figure 5 but for initial a e = &i = 1. Dimensionless time is displayed in 
panels (c,d). 
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where Mdisp 1S given by (42), M s h ear by (44), and function H(x) has the same properties as O(x) 
and We have chosen 10 



n(x) = exp 



(2- 



1.45 



(B6) 



Comparison of theoretical mass accretion rates obtained using (B5) with those derived from the 
orbit integrations in the intermediate dispersion regime is shown in Figure 7c, d. 



C. Embryo accretion rate in an inhomogeneous planetesimal disk. 

The accretion rate of any massive body immersed in a planetesimal disk depends sensitively on 
the dynamical state of the disk. Detailed studies of various accretion regimes in homogeneous disks 
were carried out by a number of authors (Greenzweig & Lissauer 1990, 1992; Dones & Tremaine 
1993) and we will use some of their results. 

We will concentrate on one very important regime: when the maximum planetesimal impact 
parameter l max necessary for accretion to take place is much smaller than the disk vertical scale- 
height. If this condition is not fulfilled the inclination of planetesimals must be unrealistically small 
(Dones & Tremaine 1993), 

R 

d-i<P, p= (CI) 
Rh 

The dimensionless parameter p is the ratio of the embryo's radius R e to its Hill radius Rh- It is 
independent of the mass of the accretor and is typically rather small in the Solar system (< 10~ 2 ). 

In the dispersion-dominated regime, when a^+af 3> 1 [this is the 3-dimensional high-dispersion 
regime of Dones & Tremaine (1993)], analytical treatment of the accretion rate is possible, because 
the two-body approximation is valid in this velocity regime (see Paper II). In this approximation 
the accretion cross-section is 



S = irRl 



l 2G(M e + m) 



R P .v 



2 



(C2) 



where M e and m are accretor and planetesimal masses, and v is the approach velocity of the 
planetesimal particles at infinity. The second term in this formula takes into account gravitational 
focussing, which is very important in planetesimal disks. 

We proceed in a manner analogous of that of Paper II (see §4). We first fix H, e, and i of 
planetesimals. Only planetesimals with \H\ < e can experience close encounters with the embryo 
and be accreted. In addition, planetesimals on horseshoe orbits cannot be accreted either, thus we 



10 This choice of n(x) certainly depends on the accuracy of formula (C9) and our specific value of H co u, see Appendix 

C. 
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put another restriction on H: \H\ > hh s , with hh s given by (22). Also, following the discussion 
in §3.3, we need to distinguish the value of the semimajor axis difference immediately before the 
encounter H' from its value at infinity, H. The conservation of planetesimal flux means that the 
surface number densities of planetesimals immediately before the encounter N' and far from the 
embryo N are related via 

N'{H')\H'\dH' = N (H)\H\dH . (C3) 

If the dependence H'(H) is given by our simple approximation (27) then dH' /dH = H/H' and 
N'(H'(H)) = N(H); of course this would not be true for any other relation between these quantities. 

The total number of planetesimals with semimajor axes between H and H + dH, eccentricities 
and inclinations in the intervals [e, e + de] and + di] correspondingly passing the accretor in a 
unit of time is 

o t->2 

5N = -n-^N(H)^j{e,i)\H\dHdedi, (C4) 

where ip(e, i) is a distribution function of eccentricities and inclination for which we assume a 
Gaussian form (3) following Paper II, and N is dimensionless. The approach velocity of accreting 
planetesimals is [equation (60) of Paper II] 

v 2 = n 2 R 2 H [e 2 + i 2 - (3/4)(H') 2 ] . (C5) 

We have used H' in this formula because calculation of scattering in the two-body approximation 
uses values of the orbital parameters immediately before the encounter. 

The accretor will absorb all planetesimals which have impact parameters at infinity I smaller 

than 



, 2G(M e + m) 

Imax = Re\j i- -\ ^""^2 > l^ 6 ) 



a condition from which (C2) follows. Then, using equation (A7) we may write that the fraction of 
5N which can get accreted is 



fa - P(H',l < Imax) = -^Z^^ 7TB~. ~ ^ 777, ( C7 ) 



2t, 

3vr R\ nR H i\H'\y/e 2 - (H>) 2 
Then the mass accretion rate due to the population of field particles with mass m is 



oo oo oo 

M = m j dH j j dediil>(e, i)f a 



d5N 

~dW 



\H\ 



oo 

?2 



MR? f ,tt f l\ T/TT , \H\ ip(e,i)dedi 



J f dH J f j N(H)d -^ 



™ 2 e J J J \H'\ ~Ue 2 - {H'f 

-oo \H\ 



(C8) 
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where we are using the following notation: v(e,i) = v/(QRh), and p is given by (CI). 

We carry out the integration over de and di in a way similar to our calculation of scattering 
coefficients in Appendix A of Paper II. As a result we obtain equations (42) and (43). Note that (42) 
automatically takes into account the transition between strong and weak gravitational focussing 
regimes which takes place at a e ,ai ~ p" 1 / 2 » 1. We have numerically compared (42) applied to a 
homogeneous planetesimal disk (and H' = H) with the corresponding analytical accretion rates of 
Greenzweig & Lissauer (1992) and found them to agree. The conversion between H' and H turns 
out to be a significant ingredient for the accuracy of the accretion rates; the use of simple H' = H 
prescription instead of (27) leads to « (10 — 20)% bigger discrepancy with our numerical results. 

In the shear-dominated regime ["intermediate dispersion, strong gravity" case of Dones & 
Tremaine (1993)] the situation is different. The approach velocity of planetesimals is always ~ 
£IRh- Then the mass accretion rate can depend on only the vertical (and not horizontal) velocity 
dispersion since a; L determines the disk thickness and local density of planetesimals. Using simple 
scaling arguments and orbit integrations to fix constant coefficients Dones Sz Tremaine (1993) have 
demonstrated that in the shear-dominated regime with the ratio of vertical to horizontal velocity 
dispersions dija e = 0.5 the accretion rate in homogeneous disk is given by 

M^5 NmnR f H - (C9) 

In our case this result is no longer applicable because planetesimal surface density is not the 
same in different parts of the disk. In fact, when computing the accretion rate in inhomogeneous 
disk it is more meaningful to use instantaneous surface number density of planetesimals on passing 
orbits which lead to collisions with the embryo. We will assume that planetesimals on orbits near 
the horseshoe-passing boundary end up colliding with the embryo in the cold regime [see Petit & 
Henon (1986) for more accurate locations of collision bands in the shear-dominated regime]; this 
means that we will be using the instantaneous surface density N inst {H co u) where H co u ps ±1.4 [see 
(20)] as a measure of surface density in equation (C9). Using results of Paper II (Appendix B) we 
can write that 

oo 

N ""'^=wJ N{H) ^rr p 

— oo 

All these considerations allow us to adopt formula (44) for the accretion rate in the inhomo- 
geneous disk with arbitrary ratio of vertical to horizontal planetesimal velocity dispersions. Of 
course, in the homogeneous disk with di/d e = 0.5 it reproduces (C9). 

In the intermediate velocity regime we smoothly interpolate between the accretion rates rep- 
resented by formulae (44) and (42), see Appendix B. To compute the accretion rate in a disk with 
a distribution of planetesimal masses one simply needs to integrate our single-mass formulae over 
the whole planetesimal mass spectrum. 
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